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Abstract.  We  derive  a  Godunov- type  numerical  flux  for  the  class  of  strictly  convex,  homoge¬ 
neous  Hamiltonians  that  includes  H ( p ,  q)  =  > J ap 2  +  bq 2  —  2 cpq,  c2  <  ab.  We  combine  our  Godunov 
numerical  fluxes  with  simple  Gauss— Seidel-type  iterations  for  solving  the  corresponding  Hamilton- 
Jacobi  (HJ)  equations.  The  resulting  algorithm  is  fast  since  it  does  not  require  a  sorting  strategy 
as  found,  e.g.,  in  the  fast  marching  method.  In  addition,  it  provides  a  way  to  compute  solutions 
to  a  class  of  HJ  equations  for  which  the  conventional  fast  marching  method  is  not  applicable.  Our 
experiments  indicate  convergence  after  a  few  iterations,  even  in  rather  difficult  cases. 
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1.  Introduction.  Hamilton  Jacobi  (HJ)  equations  have  a  rich  pool  of  applica¬ 
tions,  ranging  from  those  of  optimal  control  theory  and  geometrical  optics,  to  essen¬ 
tially  any  problem  that  needs  the  (weighted)  distance  function  [14].  Examples  include 
crystal  growth,  ray  tracing,  etching,  robotic  motion  planning,  and  computer  vision. 
Solutions  of  these  types  of  equations  usually  develop  singularities  in  their  derivatives, 
and  thus  the  unique  viscosity  solution  [6]  is  sought. 

In  this  article,  we  focus  on  the  class  of  time  independent  HJ  equations  with 
Dirichlet  boundary  condition 

JJ(x,  Vu)  =  r(x),  «|r  =  0; 

if(x,  p)  are  strictly  convex  nonnegative,  and  lim^o  R(x,  Ap)  =  0.  We  explain  our 
method  using  the  following  important  model  equation: 

(1.1)  <j>y)  =  ^Jacj>l  +  &</>2  -  2 c4>x<$>y  =  r, 

where  0  :  K.2  i — 1R  is  continuous  and  a,  6,  c,  and  r  can  be  either  constants  or  scalar 
functions;  in  the  latter  case,  H  depends  also  on  x,  defined  on  R2,  satisfying  ab  >  c2, 
a,b,r  >  0.  With  a  =  b  =  1  and  c  =  0,  we  have  the  standard  eikonal  equation  for 
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which  many  numerical  methods  have  been  developed.  This  equation  has  the  essential 
features  of  HJ  equations  with  convex  Hamiltonians,  so  that  we  can  easily  explain  our 
algorithm,  and  is  general  enough  that  fast  marching  is  not  applicable. 

In  the  following  subsections,  we  will  review  some  of  the  solution  methods  for 
the  eikonal  equation  since  it  forms  the  motivation  of  our  work.  We  then  present  a 
fast  Gauss-Seidel-type  iteration  method  for  (1.1)  which  utilizes  a  monotone  upwind 
Godunov  flux  for  the  Hamiltonian.  We  show  numerically  that  this  algorithm  can  be 
applied  directly  to  equations  of  the  above  type  with  variable  coefficients. 

1.1.  Solving  eikonal  equations.  In  geometrical  optics  [10],  the  eikonal  equa¬ 
tion 

(1.2)  \J<t>l  +  4>l  =  r{x,  y) 

is  derived  from  the  leading  term  in  an  asymptotic  expansion 

OO 

j= o 


of  the  wave  equation 


wtt  -  c2(x,  y){wxx  +wvv)  =  0, 

where  r(x,  y)  =  1/| c(x,  y)\  is  the  function  of  slowness.  The  level  sets  of  the  solution  </> 
can  thus  be  interpreted  as  the  first  arrival  time  of  the  wave  front  that  is  initially  T. 
It  can  also  be  interpreted  as  the  “distance”  function  to  T. 

We  first  restrict  our  attention  for  now  to  the  case  in  which  r  =  1.  Let  T  be  a 
closed  subset  of  R2.  It  can  be  shown  easily  that  the  distance  function  defined  by 

g?(x)  =  dist  (x,  T)  :=  min  ||x  —  p\\,  x  =  (x,  y)  €  R2, 
per 

is  the  viscosity  solution  to  (1.2)  with  the  boundary  condition 

<t>(x,y)  =  0  for  (x,  y)  gT. 

Rouy  and  Tourin  [20]  proved  the  convergence  to  the  viscosity  solution  of  an  itera¬ 
tive  method  solving  (1.2)  with  the  Godunov  Hamiltonian  approximating  ||V<^||.  The 
Godunov  Hamiltonian  function  can  be  written  in  the  following  form: 

(1.3)  HG(p_,p+,q_,q+)  =  y//max{p+,p”}2  +  max{q+ ,  q~}2 , 

where  p±  =  D±ipij,  q±  =  D±(/>ij,  D±cf>i j  =  ±((f>i± —  (f>ij)/h  and  accordingly  for 
D±<j>i,j,  and  x+  =  max(x,0),  x~  =  —  min(a:,0). 

Osher  [13]  provided  a  link  to  time  dependent  eikonal  equations  by  proving  that 
the  f-level  set  of  c/)(x,  y)  is  the  zero  level  set  of  the  viscosity  solution  of  the  evolution 
equation  at  time  t, 


fpt  =  |  |V-0|  |  =  0, 

with  appropriate  initial  conditions.  In  fact,  the  same  is  true  for  a  very  general  class 
of  HJ  equations  (see  [13]).  As  a  consequence,  one  can  try  to  solve  the  time  dependent 
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equation  by  the  level  set  formulation  [17]  with  high  order  approximations  on  the  par¬ 
tial  derivatives  [9,  18].  Crandall  and  Lions  proved  that  the  discrete  solution  obtained 
with  a  consistent,  monotone  numerical  Hamiltonian  converges  to  the  desired  viscosity 
solution  [5]. 

Tsitsiklis  [25]  combined  heap  sort  with  a  variant  of  the  classical  Dijkstra  algorithm 
to  solve  the  steady  state  equation  of  the  more  general  problem 

IIV0II  =  r(x). 

This  was  later  rederived  in  [23]  and  also  reported  in  [8].  It  has  become  known  as  the 
fast  marching  method,  whose  complexity  is  0(N\og(N)),  where  N  is  the  number  of 
grid  points.  Osher  and  Helmsen  [15]  have  extended  the  fast  marching-type  method  to 
somewhat  more  general  HJ  equations.  We  will  comment  on  this  in  a  following  section. 

1.2.  Anisotropic  eikonal  equation.  We  return  to  the  Hamiltonian  in  ques¬ 
tion:  H(p,q)  =\J ap2  +  bq2  —  2 cpq.  Writing  the  quadratic  form  as 

ap2  +  bq2  —  2cpq  =  (  p  q  )  ^  ^  j  ^  , 

it  is  easy  to  see  that  we  can  diagonalize  the  symmetric  matrix  in  the  middle  of  the 
equation  for  our  previously  noted  choices  of  a,  b,  c  and  find  a  coordinate  system  £-77 
such  that,  after  rescaling,  the  Hamiltonian  becomes 

H(p,q)  =  s/p2  +  q2. 

The  eigensystem  of  the  above  matrix  defines  the  anisotropy.  Indeed,  the  authors  in 
[15]  proposed  to  solve  the  constant  coefficient  equation  (1.1)  by  first  transforming  it 
to  (1.2)  in  the  £-77  coordinate  system. 

This  anisotropy  occurs  in  fields  such  as  ray  tracing  in  special  media,  e.g.,  crystals, 
in  which  there  are  “preferred”  directions.  Furthermore,  we  will  see  that  it  can  be  a 
result  of  considering  the  geodesic  distance  function  on  a  manifold  M  that  is  defined 
as  the  graph  of  a  smooth  function  /. 

Let  4>  be  the  distance  function  such  that 


4>{x,y) 


min 

7C  M 


ds 


and  7  connects  the  point  (x,  y)  with  the  set  T  C  M.  The  minimizing  curve  is  called 
the  geodesic,  and  </>  the  distance  function  to  T  on  M.  Moreover,  <j>  solves 


(1.4)  ||Pvv>V</>||2  =  1,  <j>\r  =  0, 


where  y,  z )  =  f(x,  y)  —  z ,  and  the  projection  operator  [4] 


Pv  v>  =  I  — 


Vt/j  0  Vt/> 

IIVV’II2 


which  projects  a  vector  onto  a  plane  whose  normal  is  parallel  to  V0.  Using  the  fact 
that  Pvj/>  is  a  projection  operator,  a  simple  calculation  shows  that 


(1.5) 

||-Pvv>V</>||2 


1  - 


fl  +  ft  +  1 


1  - 


f 

Jv 


fi  +  fi  + 1 


_  o  fxfy  ^  , 

^  f2+f2  + 


V 
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This  is  clearly  of  the  form  of  Hamiltonians  that  we  are  interested  in.  We  will  apply 
our  algorithm  to  compute  the  geodesic  distance  later  in  this  paper. 

There  are  other  approaches  that  are  designed  to  compute  distances  on  manifolds. 
For  example,  [11]  provided  an  algorithm  to  compute  the  geodesic  distance  on  trian¬ 
gulated  manifolds.  Barth  [2]  uses  the  discontinuous  Galerkin  method  to  find  distance 
on  graphs  of  functions  that  are  represented  by  spline  functions.  In  [4],  the  authors 
embed  the  manifold  as  the  zero  level  set  of  a  Lipschitz  continuous  function  and  solve 
the  corresponding  time  dependent  eikonal  equation  (1.4)  in  the  embedding  space.  As 
we  have  mentioned  in  the  previous  subsection,  the  zero  level  set  of  the  time  dependent 
eikonal  equation  at  time  t\  is  the  fy-level  set  of  the  solution  to  the  stationary  eikonal 
equations  (see  [13]).  In  [12],  the  authors  adopted  the  standard  fast  marching  method 
to  solve  the  isotropic  eikonal  equation  in  a  thin  band  of  thickness  e,  which  encloses 
the  manifold  M,  and  proved  that  the  restriction  of  the  solution  to  M  converges  to  the 
geodesic  distance  as  e  goes  to  0.  In  [21,  22],  the  authors  provide  an  ordered  upwind 
method  to  solve  a  general  class  of  static  HJ  equations.  We  will  comment  on  their 
method  in  a  later  subsection. 

1.3.  Osher’s  fast  marching  criteria.  Since  the  fast  marching  method  is  by 
now  well  known,  we  will  not  give  much  detail  on  its  implementation  in  this  paper.  In 
general,  this  involves  a  sorting  procedure  and  the  solution  of 

(1-6)  HG(D^<l>ij,I%<l>ij,Dv_<l>ij,Dv+<l>iij)  =  1 

for  (f>ij  in  terms  of  its  four  neighboring  values.  More  precisely,  the  heap  sort  strategy 
of  the  fast  marching  method  requires  a  monotone  update  sequence.  The  updated 
value  of  a  grid  node  has  to  be  greater  than  or  equal  to  those  of  the  grid  nodes  used 
to  form  the  finite  difference  stencil.  This  amounts  to  the  condition 

pHp  +  qHq  >  0, 

which  dictates  that  the  solution  be  nondecreasing  along  the  characteristics.  However, 
if  we  use  one-sided  upwind  finite  difference  approximations  for  partial  derivatives  of 
(j)  on  a  Cartesian  grid,  it  is  equivalent  to  demanding  that  the  partial  derivatives  of 
(j)  (i.e. ,  p  and  q)  and  their  corresponding  components  of  the  characteristics  directions 
(i.e.,  dx/dt  and  dy/dt )  have  the  same  sign.  Since  dx/dt  =  Hp  and  dy/dt  =  Hq,  we 
have  the  stricter  Osher’s  fast  marching  criterion 

(1.7)  pHp  >  0,  qHq  >  0. 

It  does  not  matter  whether  the  Hamiltonian  is  convex  or  not;  as  long  as  criterion  (1.7) 
is  satisfied,  a  simple  fast  marching  algorithm  can  be  applied.  But  if  the  criterion  is 
not  satisfied,  fast  marching  cannot  be  applied  to  the  problem  on  a  Cartesian  grid.  Of 
course  there  are  Hamiltonians  that  do  not  satisfy  (1.7).  In  the  class  of  Hamiltonians 
that  we  consider,  as  long  as  c  ^  0,  it  is  likely  that  the  values  of  p  and  q  differ  to  the 
extent  that  the  above  criterion  is  no  longer  satisfied.  In  light  of  criterion  (1.7),  we 
have  also  tried  to  find  directions  £(x,  y)  and  y{x,  y)  locally  in  which  pHp  >  0,  qHq  >  0. 
However,  if  one  insists  on  using  Cartesian  grids,  the  implementation  of  this  approach 
might  be  a  bit  hairy.  We  are  interested  especially  in  tackling,  over  a  Cartesian  grid, 
problems  where  the  solution  is  nondecreasing  along  characteristics  but  where  Osher’s 
fast  marching  criterion  is  not  satisfied. 


SWEEPING  METHOD  FOR  HAMILTON 


677 


1.4.  The  sweeping  idea.  Danielsson  [7]  proposed  an  algorithm  to  compute 
Euclidean  distance  to  a  subset  of  grid  points  on  a  two  dimensional  grid  by  visiting 
each  grid  node  in  some  predefined  order.  In  [3],  Boue  and  Dupuis  suggest  a  similar 
“sweeping”  approach  to  solve  the  steady  state  equation  which,  by  experience,  results 
in  an  O(N)  algorithm  for  the  problem  at  hand.  This  “sweeping”  approach  has  recently 
been  used  in  [24]  and  [27]  to  compute  the  distance  function  to  an  arbitrary  data  set 
in  computer  vision.  In  [26],  the  author  provides  some  theoretical  evidence  indicating 
that  sweeping  converges  to  an  approximate  Euclidean  distance  function,  i.e. ,  to  an 
approximate  viscosity  solution  of  |V</|  =  1  in  2d  predetermined  iterations.  We  will 
talk  about  these  iterations  in  a  later  section.  Using  this  “sweeping”  approach,  the 
complexity  of  the  algorithms  drops  from  O(NlogN)  in  fast  marching  to  O(N),  and 
the  implementation  of  the  algorithms  becomes  a  bit  easier  than  the  fast  marching 
method  that  requires  heap  sort. 

This  sweeping  idea  is  best  illustrated  by  solving  the  eikonal  equation  in  [0, 1]: 
\ux\  =  1,  u(  0)  =  u(l)  =  0. 


Let  Ui  =  u(Xi)  denote  the  grid  values  associated  with  the  uniform  grid  composed  of 
the  gridpoints  0  =  Xo  <  x\  <  ■  ■  ■  <  xn  =  1.  We  then  solve  the  discretized  nonlinear 
system 

(1.8)  \J max(max(D_Ui,  0)2,  min(H+Mj,  0)2)  =  1,  uq  =  un  =  0, 


by  our  sweeping  approach.  We  initially  set  =oo,  i  =  1, . . . ,  n  —  1.  In  practice,  oo 
can  be  replaced  by  some  number  K,  which  is  larger  than  maxje[0,i]  u.  Let  us  begin 
by  sweeping  from  0  to  1;  i.e.,  we  update  iq  from  i  =  1  increasing  to  *  =  n  —  1.  This  is 
“equivalent”  to  following  the  characteristics  emanating  from  xq.  Let  denote  the 
grid  values  after  this  sweep.  We  then  have 

u(i)  =  /  i  if  i  =  1, .  ■ . ,  n  —  2, 

*'  l  MH  =  n-l. 

Notice  that  at  i  =  n—  1,  we  actually  use  the  upwind  information  from  the  neighboring 
right  boundary  point.  Furthermore,  notice  that  already  has  the  correct  desired 
values  for  i  <  n/2  since  the  sweep  goes  from  left  to  right,  the  desired  upwind  direction 
for  these  i.  In  the  second  sweep,  we  update  iq  from  i  =  n  —  1  decreasing  to  1,  starting 
with  During  this  sweep,  we  follow  the  characteristics  emanating  from  xn.  The 
use  of  (1.8)  is  essential,  since  it  determines  what  happens  when  two  characteristics 
cross  each  other.  It  is  then  not  hard  to  see  that,  after  the  second  sweep, 

f  n  if  *  —  f  > 

Ul  1  otherwise. 

v  n 


Notice  that  the  correct  values  at  i  <  n/2  derived  after  the  first  sweep  are  un¬ 
changed,  and  new  and  correct  values  for  i  >  n/2  are  created.  In  summary,  this  simple 
iterative  algorithm  can  be  described  as  follows:  at  the  fcth  iteration,  solve 


max 


|max 


( k ) 

U )  —  U, 


(k-1) 


Ax 


,  min 


..(fc-U 

li+l 


—  U. 


(k) 


Ax 


=  1 


for  u^k:>  for  each  i  going  from  1  to  n  —  1  in  the  first  iteration  (k  =  1),  and  from  n  —  1 
to  1  for  the  second  iteration  ( k  =  2).  However,  for  more  complicated  equations  and 
boundary  conditions,  it  is  not  so  easy  to  write  down  the  equivalent  explicit  solution. 
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In  this  paper,  we  will  extend  this  sweeping  approach  to  a  class  of  HJ  equations 
that  cannot  be  solved  by  the  fast  marching  algorithm,  by  first  deriving  a  Godunov 
Hamiltonian. 

In  [21,  22],  the  authors  proposed  a  one-pass  method  that  is  based  on  a  control- 
theoretic  viewpoint.  In  principle,  they  solve  the  HJB  equation 

(1.9)  max  Vtt  •  a/(a,  x)  =  1, 

a 

where  p  =  (p,  q)  and  the  function  /(a,  x)  is  the  speed  of  motion.  This  formula  is  the 
second  Legendre  transform  taken  on  the  sphere;  see,  e.g.,  [16,  19]. 

The  idea  is  still  to  follow  the  characteristics  and  update  the  grid  value  in  a  mono¬ 
tone  sequence.  In  a  notation  similar  to  the  two  dimensional  setting  of  [21,  22],  we 
let  uQ  be  the  grid  value  we  are  updating.  To  update  ua,  we  have  to  look  for  two 
other  grid  values  ur  and  ua ,  which  are  not  necessarily  the  immediate  grid  neighbors 
of  u0-  For  example,  if  u0  is  the  grid  value  itjj,  the  immediate  neighbors  of  u0  are 
then  and  Uij- 1.  As  we  indicated  in  the  previous  subsection,  it  is 

possible  that  u0  is  less  than  all  its  immediate  neighboring  values.  We  then  need  to 
find  two  other  grid  values,  here  denoted  as  ua  and  ur,  to  form  an  upwinding  stencil. 
Then  u0  is  found  by  minimizing  a  nonlinear  expression  derived  from  (1.9),  using  the 
values  of  ur,  us,  and  /.  The  heap  sort  data  structures  are  used  in  order  to  find  ur  and 
us;  therefore,  the  complexity  is  N  log  N,  where  N  is  the  total  number  of  grid  points. 
Also,  since  ur  and  us  may  not  he  on  the  immediate  neighbors,  this  algorithm  may 
need  a  larger  region  around  the  initial  wave  front  to  get  started. 

As  one  will  see  in  the  following  section,  our  proposed  method  is  also  based  on 
following  the  characteristics.  To  update  u0,  our  method  uses  only  the  immediate 
neighboring  grid  values  and  does  not  need  the  heap  sort  data  structure.  More  impor¬ 
tantly,  our  algorithm  follows  the  characteristics  with  certain  directions  simultaneously, 
in  a  parallel  way,  instead  of  a  sequential  way  as  in  the  fast  marching  method.  The 
Godunov  flux  is  essential  in  our  algorithm,  since  it  determines  what  neighboring  grid 
values  should  be  used  to  update  u  on  a  given  grid  node  o.  At  least  in  the  examples 
presented,  we  need  only  to  solve  a  simple  quadratic  equation  and  run  some  simple 
tests  to  determine  the  value  to  be  updated.  This  simple  procedure  is  performed  in 
each  sweep,  and  the  solution  is  obtained  after  a  few  sweeps.  Our  code  is  not  much 
more  than  what  is  presented  in  section  3.2.  We  also  point  out  the  ease  of  implementing 
our  proposed  algorithm  and  its  extension  to  more  dimensions;  this  will  be  described 
in  a  sequel  paper. 

2.  A  Godunov  flux  for  strictly  convex  Hamiltonians.  By  solving  the  Rie- 
mann  problem  for  HJ  equations  (Godunov’s  procedure),  Bardi  and  Osher  [1]  proved 
rigorously  that 

(2.1)  HG(p-,p+;q-,q+)  =  ext  ext  H(p,q), 

pel\p-,p+]  qel[q-,q+ } 

where 


ext 

=  min 

p£l[a,b] 

p€[a,b\ 

ext 

=  max 

p£l[a,b] 

pe[b,a] 

HG(p-,p+;q-,q+), 
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and  I[a,  b]  denotes  the  closed  interval  bounded  by  a  and  b.  This  is  a  monotone  upwind 
flux  function,  which  implies  convergence.  Godunov’s  scheme  (1.3)  for  the  eikonal 
equation  V cf>2  +  (j)2  =  1  can  be  derived  from  the  above  formula.  It  is  one  of  the 
central  topics  of  this  paper  to  derive  an  explicit  formula  for  the  class  of  strictly  convex 
Hamiltonians  in  question.  Especially,  we  will  demonstrate  our  numerical  methods  on 
H  =  fa  a(f2x  +  b(j)y  -  2c<j>x(/>y,  c2  <  ab. 

Note  that,  in  general,  if  we  reverse  the  order  on  p  and  q  in  our  ext-ext  decision, 
the  result  might  be  different,  although  they  both  give  convergent  monotone  methods. 
However,  in  the  convex  Hamiltonian  at  hand,  the  results  are  order  independent. 

For  convenience,  we  will  also  use  Haifa, j,  fa±i,j>  fa,j± i)  t°  denote  the  evaluation 
of  our  Godunov  Hamiltonian  Hg{D^_  faj,  Di  faj,  D^_faj,  Dv,faj). 

2.1.  Derivation  of  the  flux.  In  order  to  derive  a  compact  expression  that 
satisfies  (2.1),  we  need  to  study  the  extremum  of  the  Hamiltonian  on  Ip  x  Iq  C  R2, 
where  Ip  is  a  shorthand  for  I\p_,p+\. 

The  extremum  may  occur  on  either  the  critical  points  of  H  or  the  boundary  of 
Ip  x  Iq.  Let  us  first  look  at  the  partial  derivatives  of  H,  i.e.,  Hp  and  Hq,  and  their 
zeros.  Fix  a  go;  the  extremum  of  H(p,  go)  occurs  at  either  the  critical  point  of  H{p,  go) 
(i.e.,  where  Hp  =  0)  or  the  boundary  of  I\p-,p+\.  We  denote  the  critical  point  by 
pa(q o)-  Similarly,  given  p0l  we  obtain  the  critical  point  qfapo)-  For  convenience,  we 
shall  denote  pa(qo)  by  pa  when  go  can  be  determined  from  the  context,  and  (pa,  qa)  is 
the  critical  point  of  H  such  that  Hp{pa ,  qa)  =  Hq{pa ,  qa)  =  0.  Therefore,  we  consider 
separately  H(pa,qa),  H{p_,  qa{p-)),  H{p+,qa{p+)),  H(pCT(g_),  g_),  H{pa(q+),q+), 
and  H(p±,q±)  as  possible  evaluations  of  (2.1). 

For  fixed  p ,  we  have 

(2.2)  HG{p,q-,q+)  =  H(p,  sgn  max{  (g_  -  gCT)+,(g+  -  qa)~]  +qa), 

where 


sgn max(a;,  y)  =  x+  if  max{a :+,y  }  =  x+, 
sgnmax(x,  y)  =  —y~  if  ma x{x+,y~}  =  y~ . 

The  expression  for  fixed  q  is  a  direct  analogy  to  (2.2).  It  is  easy  to  see  that  Hq{-,  •;  q~,  q+) 
is  increasing  in  g_  and  decreasing  in  q+.  By  symmetry,  Hg(p~,p+;  •,  •)  is  increasing 
in  p _  and  decreasing  in  p+ . 

Details  of  the  derivation  of  the  above  expression  are  provided  in  the  appendix. 

The  following  proposition  will  be  of  use  in  analyzing  this  introduced  Godunov 

flux. 

Proposition  1.  If  Hpp  >  0,  Hqq  >  0,  pHp  >  0,  qHq  >  0,  and  pa{ 0)  =  gCT(0)  =  0, 
then  p<j{q)  =  0  Vg,  and  qa{p)  =  0  \/p. 

Proof.  pa{q ),  by  definition,  is  the  zero  of  Hp(pa(q),q)  =  0.  We  will  write  pa  in 
place  of  pc ,  (g)  for  brevity.  This  proposition  is  then  proved  by  simple  manipulation  of 
the  definitions: 

^paH{pa,q)  =  p'„Hp(pa,  q)  +  pa{Hpp{pa,  q)p'a  +  Hpq(pa,  q)) 

d 

PcrPaHppiPm  d)  T  Hp  [Pf,  .  g) 

=  PvPa  Hpp  (Pa ,  g) 

=  0. 
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The  hypothesis  Hpp  >  0  implies  that 

pa{q)  =  0  Vg  or  p'a{q)  =  0  Vg. 

Again,  by  the  hypothesis  that  pa( 0)  =  gCT(0)  =  0,  we  can  conclude  that  pCT(g)  =  0  Vg. 
Similarly,  qa(p)  =  0  \/p.  □ 

Notice  that  if  the  Hamiltonian  is  \J p1  +  q2,  our  upwinding  expression  in  (2.2)  is 
identical  to  the  conventional  expression  max(gA,P+)-  (In  this  case,  the  sign  of  the 
second  argument  does  not  matter  since  we  are  really  evaluating  its  square  product 
in  the  eikonal  equation.)  In  fact,  we  have  the  following  corollary,  which  is  a  direct 
consequence  of  Proposition  1. 

Corollary  1.  If  Hpp  >  0,  Hqq  >  0,  pHp  >  0,qHq  >  0,  pa( 0)  =  g^O)  =  0,  and 
H{p,q)  =  H(\p\,  |g|),  then  the  Godunov  flux  can  be  simplified  to 

HG(P-,P+’,q-,q+)  =  H(max{p±,p+},max{q±,q+}). 

3.  The  sweeping  algorithms.  We  will  use  the  model  equation  (1.1)  as  a  con¬ 
crete  example  for  the  exposition  of  our  algorithm.  We  stress  here  again  that  the 
scheme  described  below  is  valid  for  a  general  class  of  convex,  homogeneous  HJ  equa¬ 
tions. 

From  the  assumption  that  the  solution  is  nondecreasing  along  the  characteristics, 

i.e., 


pHp  +  qHq  >  0, 

we  can  easily  deduce  that  the  solution  is  nondecreasing  at  least  in  either  the  x-  or 
y-  direction;  i.e.,  either  pHp  >  0  or  qHq  >  0.  Since  we  approximate  the  derivatives 
4>x{xi,j)  by  finite  differencing  using  the  neighbors  of  faj,  the  above  monotonicity 
property  translates  to  the  following  requirement  in  the  solution 

Definition  1.  Let  (fij  be  the  solution  of  Hq((/),  (f>i±ij,  4>i,j±j)  =  rjj.  We  say 
that  (f  satisfies  the  monotonicity  requirement  if 

A  TOrni{<j>i±l,j ,  <pi,j±l} ■ 

3.1.  Derivation  of  the  scheme.  Without  loss  of  generality,  we  assume  that 
r(x,y)  =  1.  Let  us  reexamine  the  equation  to  be  solved: 

(3.1)  H(p,q)  =  1, 

where 


H  :  R  x  E  ->  R. 

Equation  (3.1)  dictates  a  level  set  relation;  namely,  the  solution  is  the  1-level  set  of 
H  in  the  p-q  plane  (denoted  here  as  A).  Correspondingly,  the  solutions  of  the  HJ 
equation  with  the  Godunov  Hamiltonian 

(3.2)  HG{p+,P-\q+,q-)  =  ext  ext  H(p,q)=r 

pel\p-,p+]  q&I[q-,q+] 

satisfy  the  following  two  properties: 

•  they  are  the  intersections  of  A  and  the  set  I\p_,p+\  x  I[q_,q+\; 

•  they  are  either  the  critical  points  of  H  or  the  boundary  points  of  the  set 

i\p-,p+\  x  i[q-,q+)- 
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Fig.  1.  The  1-level  set  of  H  and  the  box  I[p—  ,p  +  ]xl[q— ,  q  +  ]. 


Figure  1  demonstrates  two  possible  configurations  of  the  intervals.  So  what  our  al¬ 
gorithm  should  do  is  find  a  suitable  value  of  <f>  on  each  grid  node  so  that  the  divided 
forward  and  backward  differences  of  <fi  at  that  grid  node  satisfy  (3.2). 

Suppose  we  are  on  the  grid  node  (i.  j).  and  it  is  determined  that 

Hc(p+,p-;q+,q-)  =  H(p_,q+)  =  l. 


Correspondingly,  for  our  model  equation  (1.1)  we  have  to  solve  the  following  quadratic 
equation: 


(3.3) 


Di,j  fyi — 1  ,j 

Ax 


1  Yi,j 

Ay 


-  2c 


Ax 


®i,j+ 1  rid 

Ay 


=  1. 


The  solution  not  only  has  to  satisfy  the  above  equation,  but  ultimately  has  to 
be  a  solution  to  (3.2),  given  its  four  neighbors  <j>i_ ij,  <j>i+ ij,  is  and  4>i,j+ 1-  The 
subfigure  on  the  right  in  Figure  1  shows  one  such  possible  configuration;  i.e. , 


Aai 


*-i  J  ^  91+1, j 


Ax 


and 


Ay 


i,j- 1  ^  YiJ+1 


Ay 


such  that 


ext  ext  H(p,q)=  min  min  H(p,q)  =  1. 

pel\p- ,p+]  q£l[q- ,q+\  p£l[p- ,p+]  q£l[q- ,q+] 

One  can,  of  course,  implement  a  tree  of  all  the  probable  cases  from  the  complete 
listing  of  that  of  the  Godunov  Hamiltonian  (2.1).  However,  we  have  a  more  straight¬ 
forward  approach  that  utilizes  the  compact  expressions  for  the  Godunov  Hamiltonian 
(2.2)  that  we  obtained  from  the  previous  section. 

Instead,  we  solve  the  equation  with  the  following  reduced  formulas  for  the  original 
Godunov  Hamiltonian: 

(3.4)  HG{p+1P--,q+,q-)  =  ext  .  fl(p  r,q), 

q€l[q-,q+\ 


-)  =  ext  H(p+,q), 
qel[q-,q+] 


(3.5) 


HG(p+,p~-,q+,q 
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(3.6)  HG(p+:p_-,q+:q_)  =  ext  H(p,q_), 

p€l\p-,p+] 

(3.7)  HG(p+,p_;q+1q_)  =  ext  H(p,q+ ), 

P€l[p-,p+ ] 

(3.8)  HG(p+,p_;q+,q_)  =  H(pa,qa). 

For  example,  in  the  first  case,  the  flux  is  equivalent  to 

H (p~ ,  sgn max{ (g_  -  qa)+ ,  (q+  -  gCT)~}  +  q„)  =  1. 

The  possible  evaluations  of  sgnmax{(g_  —  qa)+ ,  ( q+  —  q<y)~}  +  Qa  are  g_,  q+,  qa(p-), 
and  0.  We  thus  end  up  solving  the  HJ  equation  with  all  possible  arguments  for  the 
Hamiltonian. 

Suppose  we  algebraically  solve  H(p_,q+)  =  1  for  faj  and  call  the  solution  (j>can . 
We  then  compute  the  divided  differences  p±  and  q±  using  this  (f>can  in  place  of 
We  call  <pcan  valid  if  both 

H (p- ,  sgn  max{ (g_  -  qa)+ ,  ( q+  -  q<j)~}  +  q< t)  =  1, 


H (sgn max{ (p_  -pa)+,(p+  -  pa)  }  +  p0  ,q+)  =  1, 

and  4>can  satisfies  the  monotonicity  requirement  (Definition  1). 

Finally,  we  set  to  be  the  minimum  of  those  in  the  set  of  all  valid  candidate 
solutions  4>can  obtained  from  using  all  the  possible  combinations  of  the  arguments  of 
H.  This  is  motivated  by  the  first  arrival  time  interpretation  of  the  function  </>. 

In  essence,  we  are  solving  for  the  central  value  in  the  Godunov  Hamiltonian  in 
terms  of  its  four  neighbors.  It  is  well  known  and  easy  to  show  that  any  monotone 
Hamiltonian,  let  alone  Godunov’s,  is  a  monotone  function  of  this  value.  For  these 
Hamiltonians,  this  value  goes  from  — oo  to  +oo.  Thus  there  is  always  a  unique  solution. 

Definition  2  (sweeping  iteration).  A  compact  way  of  writing  this  sweeping 
iterations  in  C/C++  is  the  following: 
for (sl=-l ; sl<=l ; sl+=2) 
for (s2=-l ; s2<=l ; s2+=2) 

f or (i=(sl<0?nx: 0) ; (sl<0?i>=0 : i<=nx) ; i+=sl) 
for(j=(s2<0?ny : 0) ; (s2<0?j>=0: j<=ny) ; j+=s2) 
update  f>ij. 

3.2.  The  algorithm.  For  the  brevity  of  the  algorithm,  we  define  respectively 
hci(p,q-,q+)  :=  sgn max{ (g_  -  <?CT(j?))+,  (<?+  -  q*(p))~}  +  q<r{p), 

hG2{p-,P+,q)  ■=  sgn max{ (p_  -  P*(q))+ ,  (p+  ~  P*(q))~}  +  Pa(q), 
where  qa(p)  =  pc/b  and  pa{q)  =  qc/a. 

Algorithm.  We  assume  that  <j>(i,j )  is  given  the  exact  values  in  a  small  neighbor¬ 
hood  of  T.  We  denote  this  neighborhood  Nbd(T).  We  initialize  <j)  by  setting  <j>{i,j)  = 
to  oo.1  We  begin  by  computing  (j) (y*  for  n  =  1. 

Do  the  following  steps  while  \\4>^n'>  —  >  6  (6  >  0  is  the  given  tolerance): 

1  Notice  that  we  only  need  to  use  a  large  value  in  actual  implementation. 
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For  each  grid  point  (i,j)  visited  in  the  sweeping  iteration,  if  Xjj  ^  Nbd(r), 
do  the  following: 

(a)  For  (, sx,sy )  =  (-1, 1),  (-1,  -1),  (1,  -1),  and  (1, 1) 

i.  Solve 


H  \ 


&x‘  ($tmp(&x  >&y)  4)  &xij')')  Sy ‘(^imp^Sx  4*  ®2/))  J  


dx 


dy 


=  r(ij) 


for  (frtmpi  sxi  sy)’ 

ii.  Let 


and 


p(,SX,  Sy )  - 


q(sx,sy)  = 


Sx  ■  {(/> tmp(sx,Sy )  -  sx,j)) 


dx 


sy  ■  (</)tmp(Sx,Sy)  -  ^ n)(i,j  -  Sy)) 

dy 


iii.  Let  TGi (sx.  sy)  be  the  logical  evaluation  of  the  equality 

H(p(sx,sy),hG1(p(sx,sy),q(sx,  l),q(sx,-l)))  =  r(i,j), 
and  TG2{sx,  sy)  be  that  of 

H(hG2{p{l,  sy),p(- 1,  sy),  q{  $xi  Sy  )U{  Sxi  sv))  =r(i,j). 

iv.  Let  M(sx,sy)  =  (f>tmp(sx,sy)  -  min -  sx,j),<j>^{i,j  -  sy)). 

v.  If  Tgi(sx,  Sy),  TG2(sx,sy)  are  true  and  M(sx,  sy)  >  0,  add  4>tmP(sx,sy 

to  the  list  phi_candidate  . 

(b)  For  (sx,sy)  =  (1,0),  (-1,0) 

i.  Solve 


H 


(  Sx-(4’tmp(Sx,0)-4’<'n)('i-Sx,j))  Sx-(4ltmp(sx,0)-4l(-n',(i-Sx,j))  c\  _  (  ■  -\ 

y  dx  ’  dx  b) 


for  (bfmj) (sx ■  Sy) 

ii. 

Compute  p(sx,  t 

s y)  and  q(sx, 

iii. 

Evaluate  TGX  (s; 

Cl  Sy). 

iv. 

If  Tci(sx,sy)  is 

true  and  M (. 

phi_candidate 

For 

(SX  ,  Sy)  =  (0,  1) 

and  (0,  —1) 

i. 

Solve 

TT  (  Sy(cptmp(0,Sy)-4>,n)  dd-Sy))  C  S  y  ■  (  0t  rnp  (0 ,  Sy  )  -  <p ( (  i,  j  -  S  y  )  )  \  _  /•  .X 

n  (  dy  a’  dy  /  —  '  v‘>  J) 

for  <f>tmp(sx,  Sy). 

ii.  Compute  p(sx,  sy)  and  q(sx,sy),  following  the  definition. 

iii.  Evaluate  TG2(sx,  sy). 

iv.  If  TG2 (sx,  sy)  is  true  and  M(sx,  sy)  >  0,  add  <ptmp(sx,  sy)  to  the  list 

phi_candidate . 

(d)  Let  4>min  be  the  minimum  element  of  phi_candidate. 


<t>^n\i,j)  =  min(^(rl)(f,  j),<j)min)- 
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(e)  Clear  phi_candidate . 

2.  Set  n  =  n  +  1;  go  back  to  step  1. 

As  described  in  the  previous  section,  we  have  to  solve  the  HJ  equation  with  all  pos¬ 
sible  arguments  for  the  Hamiltonian  and  take  the  minimum  of  those  in  the  set  of  all 
valid  candidate  solutions.  The  possible  arguments  of  the  Hamiltonian  consist  of  the 
forward/backward  differences  of  4>  and  the  critical  points  centered  at  each  grid  node. 
In  the  above  algorithm,  this  set  of  all  possible  arguments  is  indexed  by  {— 1,0,  l}2. 
Therefore,  by  X(— 1, 1)  we  denote  the  quantity  X  that  is  computed  using  H(p_,q+). 
The  number  0  encodes  the  cases  of  critical  points.  For  example,  4>tmP{—  1, 1)  denotes 
the  roots  of  the  quadratic  equation  formed  by  H(p_,q+)  =  r;  </>tmp(  1, 0)  denotes  that 
of  H(p+1qa(p+)). 

We  remark  that  in  the  case  of  c  =  0,  our  algorithm  is  equivalent  to  what  is  used 
in  the  fast  marching  method  under  the  Rouy-Tourin  formula  (1.3).  Secondly,  in  our 
numerical  implementation,  we  put  a  threshold  value  in  the  evaluations  Tq i  and  Tq 2 
for  numerical  accuracy  reasons. 

4.  Examples.  Proposition  1  and  Corollary  1  show  the  equivalence  of  the  Go¬ 
dunov  flux  derived  in  this  paper  to  the  one  commonly  used  in  the  fast  marching 
applications.  The  use  of  this  sweeping  approach  with  the  Godunov  flux  (1.3)  has 
been  reported  in  [24,  26]  for  eikonal  equations;  we  will  not  repeat  those  examples  in 
this  paper.  Instead,  we  present  results  of  our  algorithm  applied  to  our  model  equation. 

4.1.  Quadratic  Hamiltonians  \J ap2  +  bq2  —  2 cpq,  ab  >  c2,  a,b  >  0. 

In  each  of  the  following  examples,  we  compute  the  difference  in  the  approximations 
in  each  successive  iteration,  i.e. ,  ||0n+1  —  (f>n\\ li;  and  say  that  the  iterations  have 
converged  if  this  distance  is  less  than  eAx,  where  e  >  0  and  Ax  is  the  grid  size.  In 
the  examples  presented  in  this  paper,  we  simply  set  the  threshold  to  be  10-10.  Notice 
also  that  the  set  T,  on  which  cj)  =  0,  is  either  a  rectangle,  an  L-shaped  piecewise  linear 
object,  or  a  set  of  isolated  points.  The  reader  can  identify  their  location  easily  from 
the  figures. 

We  started  out  by  testing  our  algorithm  on  constant  coefficient  cases.  In  the  case 


Fig.  2.  A  sweeping  result  after  2  sweeping  iterations  on  a  50  x  50  grid.  The  initial  boundary 
is  a  single  point  in  the  center,  a  =  1.0,  b  =  1.0,  c  =  0.9. 
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Fig.  3.  a  =  1,  b  =  1,  c  =  0,  with  a  more  oscillatory  r(x)  =  2.1  —  cos(47r 2xy),  on  a  200  X  200 
grid;  convergence  is  reached  in  7  sweeping  iterations.  The  subplot  on  the  left  is  the  contour  of  the 
solution  started  with  the  square  in  the  center.  On  the  right  is  the  graph  ofr(x).  Level  curves  with 
step  0.02  are  plotted. 


Fig.  4.  (A  very  degenerate  case)  a  =  0.375,  b  =  0.25,  c  =  0.29,  with  a  more  oscillatory 
r( x)  =  (2.1  —  cos(47r2xj/))/4.0,  on  a  100  X  100  grid.  Notice  that,  in  this  case,  ab  =  0.0938  is  barely 
greater  than  c2  =  0.0841.  The  contour  of  the  solution  is  plotted.  Convergence  is  reached  at  43 
sweeping  iterations. 


of  a  =  5,  c  =  0,  we  have  solutions  that  match  the  fast  marching  solutions.  Figure  2 
shows  a  result  of  a  computation  of  the  anisotropic  case  in  which  a  =  b  =  1,  c  =  0.9. 
This  is  our  first  example  in  which  the  fast  marching  method  is  not  applicable. 

Next  we  apply  the  sweeping  algorithm  directly  to  cases  in  which  the  coefficients 
of  the  quadratic  Hamiltonian  or  the  right-hand  sides  are  not  constant.  Figure  3 
shows  a  computational  result  on  a  constant  coefficient  isotropic  Hamiltonian  and 
rather  oscillatory  forcing  function.  The  rectangle  in  the  middle  is  the  set  T.  Figure  4 
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Fig.  5.  a  =  1,  b  =  1,  c{x,y)  =  0.9sin(5-7ra;),  and  r{x,y)  =  1,  on  a  50  X  50  grid.  Convergence 
occurs  after  10  iterations. 


Fig.  6.  a  =  1.5  +  sin(57ra;),  b  =  1,  c  =  —0.6,  on  a  50  X  50  grid.  Convergence  occurs  after  10 
iterations. 


shows  a  computational  result  for  a  very  anisotropic  case.  We  notice  that  the  num¬ 
ber  of  iterations  needed  for  convergence  seems  to  depend  on  the  anisotropy  of  the 
Hamiltonian  and  also  on  how  oscillatory  the  forcing  term  is.  Figures  5,  6,  and  7  show 
results  obtained  from  variable  coefficient  Hamiltonians  with  constant  and  variable 
forcing  function  r(x,y). 

4.2.  Examples  of  distance  on  manifolds.  We  now  apply  our  sweeping  algo¬ 
rithm  to  compute  the  geodesic  distance  on  manifolds  that  are  the  graphs  of  certain 
functions.  Given  a  function  f(x,y),  with  graph  z  =  f(x,y),  we  compute  the  coeffi¬ 
cients  a(x,y),  b(x,y),  and  c(x,y)  according  to  (1.5)  and  apply  our  algorithm  directly 
to  the  corresponding  HJ  equation.  We  first  test  the  algorithm  on  a  half-sphere  with 
radius  one.  Figures  8  and  9  show  the  equidistance  lines  to  one  and  two  seed  points, 
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Fig.  7.  a  =  1.5  +  sin(5-7r:r),  b  =  1,  c  =  —0.6,  and  r(x,  y)  =  2.1  +  cos(&-Kxy ),  on  a  100  X  100  grid. 
Convergence  occurs  after  10  iterations. 


Fig.  8.  This  is  an  example  of  the  distance  on  a  half-sphere.  The  sweeping  algorithm  was  applied 
to  the  graph  of  f(x,y)  =  yT.O  —  (x2  +  y2),  with  0(0, 0)  =  0  as  boundary  condition,  on  a  100  X  100 
grid. 


respectively.  Figures  10,  11,  and  12  show  similar  computation  results  applied  to 
somewhat  more  oscillatory  manifolds.  As  we  expected,  more  sweeping  iterations  are 
required  for  convergence. 

4.3.  Grid  effects.  We  first  perform  a  rotation  of  the  coordinate  system.  We 
represent  this  by 


( x ,  y)  {x,  y) 


and  let 


(«,6,c)  (a,  b,  c). 
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Fig.  9.  This  is  an  example  of  the  distance  on  a  half-sphere.  The  sweeping  algorithm  was  applied 
to  the  graph  of  f(x,  y)  =  y/l.O  —  ( x 2  +  y2),  with  two  seed  points.  Convergence  is  reached  after  2 
sweeping  iterations. 


Viti?  ■  x  pi 


Fig.  10.  The  distance  contour  from  the  seed  point  (0,0)  on  the  graph  of  f(x,y)  = 
cos(2-7nr)  sin(2-7ry),  on  a  100  x  100  grid.  Convergence  occurs  after  9  iterations. 


To  study  the  grid  effects  of  our  sweeping  algorithm,  we  set  u  =  0  on  a  rotated 
square  whose  sides  do  not  align  with  the  grid  lines.  Comparing  the  results,  shown  in 
Figure  13,  we  see  that  the  second  picture,  concentrating  especially  on  the  diamond¬ 
shaped  contour  in  the  middle,  indeed  shows  grid  effects  compared  to  the  first  picture. 
However,  with  further  grid  refinement,  as  shown  in  the  third  picture,  grid  effects 
become  unnoticeable,  and  the  solution  from  our  sweeping  algorithm  accurately  ap- 
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Fig.  11.  The  distance  contour  from  the  seed  point  (0,0)  and  (—0.8, —0.5)  on  the  graph  of 
f(x,y)  =  cos(27ra)  sin(27ry),  on  a  100  X  100  grid.  Convergence  occurs  after  11  iterations. 


Fig.  12.  The  distance  contour  from  the  seed  point  (0,  0)  on  the  graph  of  fix,  y)  =  cos(27r:r  — 
7r)  sin(27r;/  —  tt/2),  on  a  100  X  100  grid.  Convergence  occurs  after  9  iterations. 


proximates  the  exact  solution. 

4.4.  Comparison  with  the  time  marching  solutions.  We  use  the  first  order 
Runge-Kutta-  Lax-Friedrichs  method  [18]  to  discretize  the  following  equation  and 
march  to  steady  state: 

4>t  +  sgn((j)(x,y))(H(x,y,4>x,4>y)  -  r(x,y))  =  0, 


(4.1) 
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Fig.  13.  Anisotropic  case  with  a  point  source  at  (0,0).  a  =  1,6  =  l,c  =  0.9  and  a  =  1.70365 
6  =  0.296352,  and  c  =  —0.561141,  on  50  X  50  and  100  X  100  grids.  Convergence  occurs  after  2 
iterations. 


Table  1 

Comparison  of  the  time  marching  and  sweeping  solutions  to  the  example  shown  in  Figure  12. 


dx  =  2/50 

2/100 

2/200 

2/400 

2/800 

110- 

0IUi 

2.85423 

1.83377 

1.04008 

0.56206 

0.295738 

110- 

0l|oo 

1.03825 

0.708986 

0.436469 

0.246439 

0.133858 

where  cj>(x,y,t  =  0)  =  <p(x,y)  =  0  for  (x,y)  £  T  and  <f>  is  the  solution  obtained  from 
the  sweeping  algorithm. 

We  remark  that  solving  (4.1)  is  by  no  means  a  practical  method  for  solving  the 
steady  state  equation.  Thousands  of  iterations  are  required  for  steady  state,  even  if  we 
take  <p  as  the  initial  Cauchy  data.  We  use  it  only  to  verify  the  validity  of  our  algorithm. 
Secondly,  the  solutions  of  (4.1)  suffer  from  excessive  smearing  due  to  the  numerical 
viscosity  introduced  by  the  Lax-Friedrichs  method.  As  a  consequence,  (j>  does  not 
match  well  with  <fi  on  coarse  grids.  The  reader  can  compare  Figure  14  with  Figure  12, 
for  example.  However,  we  do  see  that  \\(f>  —  (j)\\  decreases  with  the  refinement  of  the 
grid  size;  see  Table  1  and  Figure  14.  We  remark  that  higher  order  approximation 
schemes  such  as  RK3-WEN05  will  greatly  reduce  the  numerical  viscosity;  the  reader 
is  referred  to  [18].  Our  purpose  here  is  only  to  show  that  the  sweeping  approximations 
converge  to  the  viscosity  solution. 
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Fig.  14.  Steady  state  of  the  time  marching  on  a  100  X  100  and  800  X  800  grid. 


Table  2 

A  numerical  convergence  study  of  the  sweeping  algorithm  applied  to  the  graph  of  f(x,y)  = 
—  ( x 2  +  y2),  with  0(0,0)  =  0  as  boundary  condition  on  the  domain  [—0.7,  0.7]  X  [—0.7,  0.7]. 


dx  =  1.4/200  1.4/400  1.4/800  1.4/1600 


110- 011-Ll 

0.0138803 

0.0079927 

0.00453004 

0.00253513 

rate 

0.796 

0.819 

0.84 

4.5.  Numerical  convergence.  Since  we  can  easily  compute  the  geodesic  dis¬ 
tance  on  a  sphere,  we  will  use  it  as  an  example  to  show  numerical  convergence  of  our 
algorithm.  A  distance  contour  plot  is  shown  in  Figure  8.  Table  2  shows  a  numerical 
convergence  of  order  1.  We  have  also  noticed  that  the  number  of  iterations  needed 
for  the  L\  difference  of  the  approximations  in  each  successive  iteration  to  decrease 
below  the  given  tolerance  seems  to  be  bounded  independently  of  the  grid  size.  This 
number  seems  to  depend  on  the  anisotropy  (c2/a6),  the  forcing  function  r,  and  the 
configuration  of  the  interface  T. 
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5.  Conclusion.  In  this  article,  we  studied  a  fast  method  for  solving  a  class  of 
time  independent  HJ  equations  with  Dirichlet  boundary  conditions.  The  Hamilto¬ 
nians  of  interest  are  homogeneous  and  convex.  This  fast  method  combines  the  idea 
of  tracing  the  characteristics  with  Godunov  construction  and  Gauss-Seidel  iterations 
with  smart  choices  of  different  updating  sequences.  In  particular,  we  discussed  some 
important  properties  of  the  Hamiltonian  H  =  \J ap2  +  bq 2  —  2 cpq,  c2  <  ab,  and  the 
corresponding  HJ  equations.  By  the  simple  structure  of  the  convexity,  we  derived  a 
compact  expression  for  the  Godunov  Hamiltonian  that  involves  taking  extrema  of  the 
Hamiltonian  in  relation  to  the  evaluations  of  the  derivatives  of  the  solution.  With 
our  compact  Godunov  flux,  the  complexity  of  evaluating  the  Godunov  Hamiltonian 
is  reduced  to  only  eight  cases  in  two  space  dimensions.  We  then  incorporated  the 
expression  into  a  simple  Gauss-Seidel-type  iteration  procedure.  We  have  produced 
some  computational  results  using  this  algorithm.  In  particular,  we  have  applied  our 
algorithm  to  compute  geodesic  distances  on  graphs  of  functions.  This  is  of  some  im¬ 
portance  since  people  are  interested  in  finding  the  geodesics  on  terrain- like  manifolds. 

We  also  remark  that  this  Godunov-flux  sweeping  approach  can  be  extended  to 
higher  dimensional  cases.  We  are  currently  preparing  another  paper  on  this  subject. 

Our  experience  shows  that  the  number  of  iterations  needed  depends  on  the  amount 
of  anisotropy  and  the  nature  of  the  forcing  function.  Under  normal  nondegenerate 
circumstances,  experience  shows  an  O(N)  complexity  for  convergence,  where  N  is 
the  number  of  grid  points.  Recently,  in  [26],  the  author  provided  some  theoretical 
evidence  on  the  bound  of  the  number  of  iterations  for  isotropic,  homogeneous  eikonal 
equations.  This  points  out  a  future  research  direction  of  bounding  the  number  of 
sweeping  iterations  needed  for  convergence  in  relation  to  the  anisotropy. 

6.  Appendix. 

6.1.  Derivation  of  the  flux  for  homogeneous  convex  Hamiltonians.  To 

obtain  the  formula  used  earlier  in  this  paper,  we  simply  verify  its  equivalence  to  the 
following  cases,  which  rely  only  on  the  convexity  of  H: 

P-  <  p+,  and  q-  <  q+  : 

Hq  =  min  min  H  (p,  q) . 

p6[p-,p+]  qe[q-, q+] 


•  If 

Qa 

e  [ q~,q+ ], 

- 

Pa  <P- 

< 

P+, 

H(p-,qa), 

- 

p-  <p+ 

< 

Pa, 

H{p+,qa), 

- 

P-  <Pa 

< 

P+, 

H(pa,  qa). 

•  If 

Q<7 

<  Q-, 

- 

Pa  <P- 

< 

P+, 

H  (p-  ,q~), 

- 

P-  <P+ 

< 

Pa, 

& 

+ 

tcj 

- 

P-  <  Pa 

< 

P+, 

H{pa,  <?-)• 

•  If 

Qa 

>  q+, 

- 

Pa  <P- 

< 

P+, 

H(p-,q+), 

- 

P-  <P+ 

< 

Pa, 

H(p+,q+), 

— 

P-  <  Pa 

< 

P+, 

H(pa,q+ )• 

p-  <  P+, 

,  and  q _ 

> 

9+  : 

Hg  =  min  max  H(p,q)  =  min  nra x{H(p,q-),H(p,q+)}. 

p+\p-,p+\q+[q+,q~]  P6b-,P+] 
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•  If  q<j  <  q+, 

-  pa  <P-  <p+ ,  H(p_,q_), 

~  P-  <P+  <P, 7,  H(p+,q_), 

-  P-  <Pa  <P+ ,  H(p  01 Q—)' 

•  If  qa  >  q- , 

-  pa  <p-  <p+,  H(p_,q+), 

-  P-  <P+  <P, 7,  H{p+,q+), 

-  P-  <Pa  <P+,  H{p  O')  q+)- 

•  If  q+  <  qa  <  <?_ , 

~  (q<7  -  q+)  >  (q-  -  qa),  H(-,q+), 

-  ( Qa  ~  q+)  <  (q-  -  qa),  H(-,q-). 

P-  >  and  q-  >  q+  : 


Hq  =  max  max  H  (p,  q) . 
pe[p+,p-]  q£[q+,q-] 


•  If  qa  >  q~, 

-  Pa  >p~,  H(p+,q+), 

-  Pa  <P+,  H{p_1q+). 

•  If  qa  <  q+, 

-  pa  >p~,  H(p+,  q~), 

-  Pa  <P+,  H (p- ,(?-). 

•  If  q+  <  qa  <  q - 

-  ( qa  -  q+)  >  (q-  -  qa),  H(-,q+), 

-  ( qa  -  q+)  <  (q-  -  qa),  H(-,q~)- 

p -  >  p+,  and  q-  <  q+  : 

Hq  =  max  min  H  (p,  q) . 

pG[p+,p_]  q£[q-,q+] 


•  If  qa  €  [q-,q+], 

-  Pa  >P-,  H{p+,qa), 

-  Pa  <P+,  H(p_,qa). 

•  If  qa  <  q- , 

-  Pa  >p~,  H{p+,q_), 

-  Pa  <P+,  H (p_ ,(?-). 

•  If  qa  >  q+, 

-  Pa  >p_,  H(p+,q+), 

~  Pa  <P+,  H(p_,q+). 
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